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Abstract 

A Monte Carlo program for the simulation of electromagnetic cascades initiated by high-energy photons and electrons 
interacting with extragalactic background light (EBL) is presented. Pair production and inverse Compton scattering 
on EBL photons as well as synchrotron losses and deflections of the charged component in extragalactic magnetic 
fields (EGMF) are included in the simulation. Weighted sampling of the cascade development is applied to reduce 
the number of secondary particles and to speed up computations. As final result, the simulation procedure provides 
the energy, the observation angle, and the time delay of secondary cascade particles at the present epoch. Possible 
applications are the study of TeV blazars and the influence of the EGMF on their spectra or the calculation of the 
contribution from ultrahigh energy cosmic rays or dark matter to the diffuse extragalactic gamma-ray background. As 
an illustration, we present results for deflections and time-delays relevant for the derivation of limits on the EGMF. 
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1. Introduction 

The Universe is opaque to tiie propagation of y-rays witii energies in the TeV region and above 1 1]. Such photons 
are absorbed by pair production on the extragalactic background Hght (EBL) |t2|-|5|], consisting mainly of infrared 
light and the cosmic microwave background (CMB). As a result the photon flux at energies E > 10 TeV from distant 
sources as e.g. blazars is significantly attenuated on the way from the source to the Earth. High-energy photons are 
however not really absorbed but initiate electromagnetic cascades in the intergalactic space, via the two processes 

r + 76 ^ + (1) 

+yh + y . (2) 

The cascade develops very fast until it reaches the pair creation threshold afl imin = AEySy — 4m^ with Sy as the 
characteristic energy of the background photons yt,. Electron^ continue to scatter on EBL photons in the Thomson 
regime with an interaction length of a few kpc, producing photons with average energy 
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using By - 2.7 TcMB ~ 6.3 x 10^^ eV as the typical energy of CMB photons. 

The resulting shape of the energy spectrum of the diffuse photon flux Jy can be estimated analytically |@] for a 
monochromatic background as 

( K(Ey/E,r^'^ for Ey<E,, 

JyiEy) = \ KiEylE^)-^ fOV E^<Ey< E 
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V for Ey > Ejnin ■ 

Here, iSmin - m^Jsy is the threshold energy for parr-production, while Ey < E^ is the energy region where the number 
of electrons remains constant. Since the last generation of e^e" pairs produced share the initial energy equally, 
Ee - E,^in/2, this transition energy is given by E^ = 4eyE^/(3ml) = iiinin/3. Thus for a monochromatic background 
the plateau region characterised by an 1 /Ey spectrum extends only over one third of an energy decade. 

A better analytical description of the cascade development in the EBL uses a dichromatic photon gas, with scmb = 
6.3 X 10 eV and sir = 1 eV as typical energies for the CMB and the (second) peak of the IR background, respectively. 
Below one half of the threshold energy of pair production on the IR, E^ ~ £min,iR/2 - m^/(2siR) ^ 1.3 x 10" eV, 
the number of electrons remains constant. In the intermediate regime, /imin.iR < £min,CMB, electrons are Compton 
scattering on CMB photons in the Thomson regime, while photons are still producing e^e^ pairs on IR photons. Thus 
the energy E^ below which no additional electrons are injected in the cascade is given by 

„ 4 fiCMB-Ef 1 ecMB „ 

£x = - 5 = - £min.IR ~ 50 MeV . (4) 

3 3 fiiR 

Let us now compare how well this qualitative picture agrees with the cascade spectrum calculated with our Monte 
Carlo simulation. Figure [T| shows the results obtained with ELMAG for the diffuse spectrum of secondary photons 
produced by photons injected with energy E = 10'^ eV at redshift z = 0.02 and z - 0.15. We note first that the slope 
of the photon spectrum below Ey^ and the obtained value of E^ agree very well for both distances with the prediction 
in the simple "dichromatic model," although the latter assumes an infinite number of interactions. In contrast, both 
the extension and the shape of the plateau region are less universal: The smaller the distance to the source, the less 
pronounced occurs the steepening of the photon spectrum from the £ ' Thomson slope towards the predicted l/E^ 
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Figure 1: The (normalised) diffuse plioton flux E^J(E) for two sources injecting photons with energy E = 10" eV at redshift z = 0.02 (— ) and 
z = 0. 15 ( — ), respectively. 



plateau. However, the deviation of the slope from the prediction, E versus E ^, is minor akeady for distances of 
~ 500 Mpc. 

An important application of electromagnetic cascades is the calculation of various contributions to the extragalac- 
tic diffuse gamma-ray background (EGRB). Since the Universe acts as a calorimeter for electromagnetic radiation, 
accumulating it in the MeV-TeV range, the measured EGRB limits all processes during the history of the Universe 
that inject electromagnetic energy above the pair creation threshold. Examples for such processes are photo-pion and 
P + ycMB — * p+e^+e^ pair-production of UHECR protons interacting with the cosmic microwave background (CMB) 
JtIi, the decay or annihilation of (superheavy) dark matter or of topological defects |8]. Another important application 
of electromagnetic cascades is the calculation of spectra from point sources as TeV blazars. If the spectra of such 
sources extend to sufficiently high energies, emitted photons interact with the EEL. The charged component of these 
cascades is deflected by extragalactic magnetic fields (EGMF), leading potentially to halos around point sources 
1211 . to delayed echos of flaring emission 11311 and influences the observed energy spectrum 1 14]. A detailed modelling 
of the electromagnetic cascade process is thus not only necessary to connect the observed energy spectra of TeV 
sources with their intrinsic spectra, but provides also information about EGMFs. 

The extremely small interaction length compared to typical source distances from hundreds of Mpc to Gpc means 
that a large number n of interaction steps has to be simulated using a Monte Carlo approach. The exponential growth 
of the number N - 2" of secondaries aggravates the computational load in a brute-force Monte Carlo approach. 
The Monte Carlo program presented here uses weighted sampling of the cascade development to reduce efficiently 
the number of secondary particles which are traced explicitly. For maximally weighted sampling, the number of 
secondaries stays on average constant as function of interaction steps. Synchrotron losses and deflections of the 
charged component in extragalactic magnetic fields (EGMF) are included in the simulation too. The version presented 
here is restricted to the limit of small deflections. 



2. Modelling of the cascade process 

2.1. Interaction rate of photons and electrons 

The interaction rate Ry{E, z) of photons with energy E at redshift z can be connected to the pair-production cross 
section crpaii (s) and the spectral density of background photons ny{E,z) as 

RyiE, z) = ^ J dE' ny(E',z) J t/yu (1 - yu) o-p^iAs) &(s - s^in) 



We use natural units, h = c = = I, throughout the text. 
^We call from now on electrons and positrons collectively electrons. 
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where we introduced the auxiHary function 



X^""« dE' 
ny(E',z). 



(5) 



(6) 



Here we have also assumed that the EBL, as any truly diffuse background, is isotropic. 

The cm. energy squared in a yy interaction is 5 = 2EE'{1 - ji) with fi - cos while the integration limits are 
given by the pair production threshold i^in = and s^^^iE) - 4£'£'max with fimax ~ 14 eV as the high energy cutoff 
of the EBL background. The well-known pair-production cross section o-^^i^is) is given by 



' pair 
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with cTxh - 87Ta^/(3ml) as Thomson cross section and /3 - s/l - 4m; /i. 

Electrons emit in the Thomson regime mainly soft photons, cf. Eq. ([3]). To speed up the simulation, we include 
therefore as discrete interactions only those which produce secondary photons above an arbitrary energy threshold 
iithr- The remaining soft interactions are integrated out and included as continuous energy loss. Thus we define the 
interaction rate ReiE, z) of an electron with energy E at redshift z as 

Re(E,z) = dE' riyiE' , z) J d^ (1 - /3n) a-cis,e)@(s - Synm(s)) 

I /-Sm«(£) 



2 I s — ml 
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with s - Ethi/E,/3 - - m^/E^, and s - ml+2EE'{l -/3ij). The integration limits are given by Sinin(e) = m^/(l -e) 
and 5max(£) - '"e +2£'£'max(l while the Compton scattering cross section crc(s, e) integrated above the threshold 
e is given by 



crc(s,fi) - - (Txh ymin — ; 

4 1 - Jmin 
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(9) 



Here, y^^^ - m^/s and y^ax - 1 - e are respectively the minimal and the maximal energy fractions of the secondary 
electron. 

In turn, for the electron energy loss per unit distance due to the emission of photons of energies E < E±i- one 
obtains 



dEjcs/tia , , 3 

-, {S,E) = jO-Thyn 

dx 4 
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(10) 



In the left panel of Fig.|2]we show the interaction rates Ri of electrons and photons at thepresent epoch as function 
of energy. The difference between the "best-fit" and the "lower-limit" EBL from Ref. [5] becomes visible in the 
interaction rate Ry of photons only in the energy range below 10'"^ eV. The interaction rate R^ of electrons is shown 
only for the "best-fit" EBL but for two different values of the threshold isthr used in the Compton scattering cross 
section, ^thr = 3 x 10"* eV and ^thi- = 3 x 10^ eV. Note that while ^thr = 3 x lO^'eV leads akeady below lO" eV 
to strong deviations from the Thomson scattering cross section, the resulting photon spectrum is influenced by the 
threshold mainly at energies in the MeV range and below, cf. Eq. (O. 
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Figure 2: Left: Interaction rate R at z = as function of the energy E for electrons with ZJth = 3 X 10* eV and £th = 3 X 10* eV and for photons 
with the "best-fit" and the "lower-limit" EBL from Ref. |5]. Right: Fluence contained inside the 95% confidence contour of the PSF of Ferrai-LAT 
as function of energy for EGMF with top-hat profile and filling factor / varying from / = 0. 1 to / = 0.9 with £n,ax = 20TeV. 



2.2. Interactions 

The modelling of yy and ey interactions starts from sampling the cm. energy squared s of the collision according 
to the integrands of Eqs. ^ and ([8]), respectively. Technically, the rejection method is used as in most other cases 
to choose s according to its probability distribution: The value of s is first sampled logarithmically in the interval 
[■Smin, Smax], then the choice is accepted with the probability proportional to s times the integrand of Eq. (|5]) and ([8]), or 
otherwise rejected. 

For given s, the energy fraction y of the lowest energy secondary lepton (electron or positron) in the pair production 
process is sampled according to the corresponding differential cross section 



oc — 

dy y 
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with /3 - -\Jl - Arn^l s. The other secondary lepton has then the energy fraction l-y. Similarly, the energy fraction y 
of the secondary electron in the inverse Compton process is sampled according to the differential cross section 



d<Tic5{s,y) 1 

oc — 

dy y 



1 + y^ 2yniin(3' - ymin)( 1 - y) 



Kl -Jmin)^ 



(12) 



with y„ 



■ mils. 



2.3. Stacking and weighted sampling 

The produced secondary particles are then subject to a weighted sampling procedure: A secondary particle car- 
rying the fraction y of the parent energy is discarded with the probability (1 - y^^-^V'^)^ or added with the probability 
yffsampk (Q fjjg stack. Depending on the choice of the sampling parameter (0 < ^sample ^1) either all the secondaries 
are kept in the cascade (csampie = 0) or only some representative ones are retained. In particular, one secondary per 
interaction is retained on average for the default value asampie = 1 . As compensation, each particle in the cascade 
acquires a weight w which is augmented after each interaction as w ^ w/y°'"^p'° . The particles in the stack are ordered 
according to their energies. After the interaction, the lowest energy particle is extracted from the stack and traced 
further in the cascade process. 

The optimal value of Csampie depends on the typical energy of the injected photons. If the latter is so low that the 
cascades consists on average of only few steps, reducing ^sample may be advantageous because the fluctuations are 
thereby reduced. 
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2.4. Synchrotron losses 

Synchrotron energy losses of electrons are accounted for in the continuous energy loss approximation using the 
interpolation formula Ii23[l . 

^ . (13) 

dx [1 + 4.8(1 + x) ln(l + l.lx) + 3.44;^r^]2/3 

with;t' - (p±/fne){B/Bci), where p± denotes the momentum perpendicular to the magnetic field and Bcr = 4.14xlO'^G 
the critical magnetic field. 

2.5. Angular deflection and time delay 

For the energies considered, Ey > MeV and E^ > lOGeV, secondary particles are emitted in the forward direc- 
tion. Thus the angular deflection of the cascade particles results from the deflections of electrons in the extragalactic 
magnetic field (EGMF). If the coherence scale of the EGMF is much larger than electron mean free path Ae - R^^, 
an elementary deflection angle of i-th electron in the cascade chain can be calculated assuming the field to be regular 
over the distance c/, travelled by the electron, 

with being the momentum component perpendicular to the local direction of the magnetic field. Inside the patch j 
of a chosen coherence length, the deflection angles /?, per electron path are summed up coherently, /3j - YiiPi- The de- 
flections angles per coherent magnetic field patch are then summed quadratically in the random- walk approximation 
to obtain as total deflection, i.e. the angle between the initial and final photons in the cascade. 



In the small-angle approximation and assuming spherical symmetry, the angle (3 between the initial and final photons 
in the cascade is related to the emission angle a and the observation angle i? as or = cf.Ref. [ll. 

As the energy of the cascade particles quickly degrades along the cascade chain, the largest contribution to p 
comes from the last electron in the chain. This allows us to approximate the corresponding geometry by a triangular 
configuration and to obtain as relation between fi and § 

X 

sinj? = — sin/?. (16) 

Here x refers to the distance from the source S to the point P where the final photon in the cascade branch has been 
created and L is the total distance between the source and the observer O. For small i?, we thus have 

j?=|sin/3. (17) 

The time delay Afgeo of photons with respect to the straight line propagation from the source is then 

Afgeo - x{l+ sin a/ nind) - L^2x(l - x/L) sin^yS. (18) 

We add to this geometrical time delay Afgeo the kinematical time delay Afkin due to velocity v < c of the electron, 
although the latter is usuafly negligible. 

2.6. Cosmology 

The connection between redshift z, comoving distance r and light-travel time f calculated for a flat Friedmann- 
Robertson-Walker universe with Qa - 0.7 and Q,„ - 0.3 is contained in the file redshift. 
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3. Programme structure 



The programme is distributed among the files moduleslOl . f 90, userlOl . f 90, initlOl . f 90, elmaglOl . f 90 
and auxlOl . f 90. The file moduleslOl . f 90 contains the definition of internal variables, mathematical and physical 
constants; for standard applications of the programme no changes by the user are needed. The file userlOl . f 90 
contains the input/output subroutines developed by the user for the desired task. An example file is discussed in 
Sec. 5. Data files of the used EBL backgrounds and the cosmological evolution of the universe are provided in 
the directory Tables. They are read by the subroutines init_EBL(myid), init_arrays (myid) and the func- 
tion aintlR(E,z) inside the file initl01.f90. Then the function w_EBL_density_tab(emiii,zz) tabulates the 
weighted background photon density ly defined in Eq. followed by the tabulation of the interaction rate in the the 
subroutine rate_EBL_tab (eO , zz , icq) and of the electron energy losses due to the emission of photons with energy 
below the threshold in eloss_thr_tab (eO , zz , icq) . 

We discuss now in more detail the subroutines and functions of the file elmaglOl . f 90 which constitute the core 
of the programme: 

• subroutine cascade ( icq, eOO,weightO,z_in) 

Follows the evolution of the cascade initiated by a photon (icq = 0) or an electron/positron (icq = +1) 
injected at redshift z_iii with energy eOO and weight weightO until all secondary particles have energies below 
the energy threshold ethr or reached the observer at z = 0. 

• subroutine angle_delay(the2,xx,rcmb,theta,dt) 

Determines the photon time-delay dt and the observation angle theta from the rms cascade deflection angle 
the2 and the photon emission point xx by the parent electron/positron. 

• subroutine interactionCeO ,xO ,zz, t .weight ,thel ,the2,xxc ,xx,dt , icq) 

Handles one interaction with background photons: determines the cm. energy sgam of the reaction via a call 
to sample_photon or sajnple_electron(eO,zz,sgam,ierr), the energy fraction z of secondaries via a 
call to the functions zpairCsgcun) or zics(eO,sgam), and stores then the secondaries calling the subroutine 
store_particle. 

• subroutine store_particle (eO , xO , zz , t , ze , weight , thel , the2 , xxc , xx , dt , icq) 

Decides if a produced secondary is stored using weighted sampling; if yes, it adds the secondary to the array 
event and re-orders the array according to the particle energies. 

• subrout ine get_part icle(eO,xO,zz,t, weight , the 1 , the2 , xxc , xx , dt , icq) 

Reads the secondary with the lowest energy out of the array event and reduces the particle counter j cmb by 
one. 

• subroutine sample_photon(eO,zz,sgam,ierr) and sample_electron(eO,zz, sgam, ierr) 
Determines the cms energy sgam of an interaction at redshift zz. 

• double precision function w_EBL_density (emin.zz) 
Determines the weighted background photon density ly defined in Eq. (|6]l. 

• function int_point(eO,xO,zz,icq) 

Finds the next interaction point for interaction with EBL photons. 

• function sigpair (sgam) and sigics (eO , sgsim) 

Calculate the pair production and inverse Compton cross section, respectively. 

• function zpair (sgam) and zics (eO , sgsun) 

Determine the energy distribution in pair production and inverse Compton scattering, respectively. 

• function zsigics (eO , sgaim) 

Calculates the electron energy losses per unit distance due to photon emission below the chosen threshold in 
Compton scattering. 
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Figure 3: Fluence contained inside tlie 95% confidence contour of the PSF of Fermi-LAT as function of the time-delay for i? = 10 G; left for 
individual time layers, right cumulative times. 



• function zloss(eO,zz) 

Interpolates the integrated energy loss due to emission of photons below the threshold. 

• subroutine rate_EBL(eO,zz,icq) 
Interpolates the interaction rates on EBL photons. 

• function eloss_syn(E,begmf ) 

Calculates the synchrotron losses according to Eq. (fTSl) . 

• function themf (eO,dx,begmf ) 
Determines the deflection angle in the EGMF. 



The file auxlOl . f 90 contains auxiliary functions, e.g. the random number generator psran from Ref. Il24l] . 



4. Example input and output 

The file userlOl . f 90 is an example file for the input/output subroutines which should be developed by the user 
for the desired task. We discuss now the example contained in the distribution. 

4.1. Example input 

The input variables specified in the module user_variables are: the choice of the EBL model (model), the 
number of injected particles (nmax), the jet opening angle of the source in degrees (th_j et), the sampling parameter 
a_smp, the energy threshold ethr for Compton scattering, and the maximal photon energy egmax. The last two 
parameters serve also as minimal and maximal energy in the energy spectra produced as output. In subroutine 
user_main(myid,nmax) the initial redshift z and the particle type icq of the injected particles is fixed. 



z = 0.14dO 
do nl=l,nmax 

call initial_particle (eO .weight) 

icq = 

call cascade (icq, eO .weight ,z) 
enddo 



! initial redshift 

generate initial energy 
(0 - gamma, +-1 - e+-) 
starts e/m cascade 



The subroutine initial_particle (eO , weight) chooses the energy and the weight of one initial particle in the 
energy range [emin , egmax] according to a broken power-law with exponents garni below ebreak, and gajii2 above. 

The magnetic field B is modeled as patches of uniform field-strength |B| of size Zcoh- The value of the coher- 
ence length is fixed by the parameter cohlnth in module user_variables, the field-strength perpendicular to the 
propagation direction in the function bemf (r) . 
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4.2. Example output 

The energy eO, the observation angle theta and the time delay dt of secondary cascade particles with weight 
weight of type icq reaching the observer at z = are recorded by the subroutine register (eO , theta , dt , weight , icq) 
and binned in various data arrays defined in the module user_result. All data arrays exists in two versions, e.g. 
spec(n_bin,0: 1) and spec_tot (n_bin, : 1). Using MPI iflsll . the former arrays contain the result of a single 
process, which are summed by call MPI_REDUCE into spec_tot (n_bin ,0:1), 

n_array = 2*n_biii 

call MPl_REDUCE(spec,spec_tot,n_array,MPI_DDUBLE_PRECISION,MPI_SUM,0, & \\ 
MPI_COMM_WORLD, ierr) ! sum individal arrays spec 

Finally, the subroutine user_output (n_max,n_proc) writes the data arrays with the results in the files contained 
in the subdirectory Data. 

The file spec_dif f contains the normalised (diffuse) energy spectra of photons and electrons, in the format EjeV, 
E^dNy/dE and E^dNJdE. The file spec_95 includes the energy spectra of photons inside and outside the 95% area 
of the point-spread function of Fermi-LAT, in the format E/eV, E^dNy/dE(§ < 1^95) and E^dNyldEid- < §95). An 
approximation to the point-spread function of Fermi-LAT is defined in the function thereg_en(en) . 

The right panel of Fig. |2] shows the energy spectra of photons arriving within the 95% area of the point-spread 
function of Fermi-LAT for different filling factors of the EGMF which can be chosen by the parameter f rac in the 
function bemf . Otherwise the default values contained in the distributed file userlOl .f 90 are used. The photon 
fluence is compared to H.E.S.S. data j 16] and upper limits from Fermi-LAT IitIi for the TeV blazar lES 0229+200. 

The files spec_95_t and spec_95_c contain the energy spectra of photons arriving within the 95% area of the 
point-spread function, with the time-delay binned in seven time intervals, t < lOyr, lOyr < lO^yr,. . .f > 10^ yr. 
The file spec_95_c is the cumulative version of the distribution in spec_95_t. Figure [3] shows the fluence inside 
the 95% PSF of Fermi-LAT for an injection spectrum dNyjdE oc E^^^^ at redshift z - 0.14 with maximal energy 
£„,ax = 20TeV. 

5. Possible extensions 

We discuss four possible applications of the simulation ELMAG and the required extensions to perform them, 
ordered by the complexity of the necessary changes and additions. 

EGRB from dark matter decays or annihilations. High-energy electrons and photons can be generated by decays or 
annihilations of sufficiently heavy dark matter particles. For instance, the annihilation mode XX e^e^ of the dark 
matter particle X with mass nix would correspond to the injection of two electrons with energy E^ = mx- The only 
necessary addition for the calculation of the resulting EGRB is a subroutine choosing the injection point according 
to the so-called boost factor B(z) which accounts for the redshift dependent clustering of dark matter in galaxies. 
Additionally, the desired fragmentation functions dNyjdE of the X particles should be included into the subroutine 
initial_particle in the case of photons from hadronic decay or annihilation modes. 

EGRB from UHECRs and cosmogenic neutrinos. The Greisen-Zatsepin-Kuzmin cutoff" is a steepening of the proton 
spectrum at the energy Ecza ~ (4 - 5) x 10'*^ eV, caused by photo-pion production on the CMB. An additional 
signature for the presence of extragalactic protons in the cosmic ray ffux and their interaction with CMB photons is the 
existence of ultrahigh energy cosmogenic neutrinos produced by charged pion decays |18], while the corresponding 
ffux of cosmogenic neutrinos from ultrahigh energy nuclei is suppressed. Photons and electrons from pion decay and 
p + ycMB p + e* + e^ pair-production lead to a contribution to the EGRB which can be used to limit cosmic 
rays (CR) models and ffuxes of cosmogenic neutrinos. To perform this task, ELMAG has to be coupled to a program 
performing the propagation of ultrahigh energy cosmic rays which provides secondary electrons and photons from 
CR interactions as input. As the communication between the two program parts is restricted to calls of the subroutine 
cascade(icq,eO, weight, z), such a combination should be straightforward. For an example where ELMAG was 
used in this context see Ref. iH. 
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Extension to 3-dimensional cascades. Going beyond the small-angle approximation requires to calculate the actual 
trajectory of electrons solving the Lorentz equation. Additionally, scalar quantities like x , xxc , eO , begmf , . . . have 
to be changed into three-dimensional vectors and the type one.event in the module stack has to be adjusted. 
Finally, the image of a three -dimensional cascade can be calculated using the method described in Ref. lEoll . For an 



illustration of possible applications of ELMAG to this problem see Ref. |l21 1 



Treatment of interactions in sources. Photons and electrons are generated often in sources containing dense photon 
fields, as e.g. near the cores of active galactic nuclei. In this case, electromagnetic cascades take place on non- 
thermal, anisotropic photon backgrounds inside the source before the escaping particles cascade on the EBL. In order 
to describe both cascades on the EBL and inside the source, subroutines as e.g. irLit_EBL have to be doubled, adding a 
corresponding subroutine iiiit_source for the photon field inside the source. All existing subroutines which depend 
on the chosen EBL (i.e. use EBL_f it) have to adapted. In particular, the rejection mechanism in subroutines as e.g. 
sample.photon has to adjusted. As a result, such an extension requires considerable work and thorough tests of the 



changed code. A short discussion of anisotropic photon fields is given in Ref. II22I1 . 



6. Summary 

We presented a Monte Carlo program for the simulation of electromagnetic cascades initiated by high-energy 
photons and electrons interacting with the extragalactic background light. The program uses weighted sampling 
of the cascade development and treats Thomson scattering below a chosen threshold in the continuous energy loss 
approximation in order to speed up computations. 

Possible applications are the study of TeV blazars and the influence of the EGMF on their spectra or the calcu- 
lation of the contribution from ultrahigh energy cosmic rays or dark matter to the diffuse extragalactic gamma-ray 
background. As an illustration for possible applications we presented results for deflections and time-delays relevant 
for the derivation of limits on the EGMF studying the spectra of TeV blazars. Other possible applications include 
e.g. the calculation of the contribution from ultrahigh energy cosmic rays or dark matter annihilations to the diffuse 
extragalactic gamma-ray background. 
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